Quantum magnetic properties of the SU(2N) Hubbard model in the square lattice: a 

quantum Monte Carlo study 



o 

(N 



O 

CO 



u 

c3 



I 

O 

o 



> 

CO 

oo 

o 

CN 



Zi Cai, 1 Hsiang-Hsuan Hung, 2,1 Lei Wang, 3 Yi Li, 1 and Congjun Wu 1,4 

1 Department of Physics, University of California, San Diego, CA92093 
2 Department of Electrical and Computer Engineering, University of Illinois, Urbana, Illinois 61801 
Theoretische Physik, ETH Zurich, 8093 Zurich, Switzerland 
^School of Physics and Technology, Wuhan University, Wuhan, Hubei 430072, China 

We employ the determinant projector quantum Monte-Carlo method to investigate the ground 
state magnetic properties in the Mott insulating states of the half-filled SU(4) and SU(6) Fermi- 
Hubbard model in the 2D square lattice, which is free of the sign problem. The long-range antifer- 
romagnetic Neel order is found for the SU(4) case with a small residual spin moment. The SU(6) 
case is a spin gapless quantum paramagnetic state with the absence of long-range antiferromagnetic, 
dimer, and charge flux orderings. Spin correlations exhibit power-law decay, which is consistent with 
a possible algebraic spin-correlation state with an odd number of fermions per unit cell. 
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Quantum antifcrromagnetism (AF) has been an im- 
portant topic of the two-dimensional (2D) strongly cor- 
related systems for decades. For the Hubbard model in 
the 2D square lattice, charge gap opens starting from an 
infinitesimal U. The low energy physics is described by 
the AF Heisenberg model. For the SU(2) case, quan- 
tum spin fluctuations are not strong enough to suppress 
AF long-range order [J Augmenting the symmetry 
to SUYTV) or Sp(27V) enhances quantum spin fluctua- 
tions [3|45( , which can be handled by the systematic 1 /TV- 
analysis. The SU(7V) spin operators can be formulated in 
terms of either bosonic or fermionic representations. The 
bosonic large-TV analysis finds gapped quantum param- 
agnetic states exhibiting various crystalline orderings [6j , 
while the fermionic one gives rise to gapless flux-type 
spin liquid states [H, @]. However, its stability remains an 
open issue. On the other hand, short-range resonating- 
valence-bond type gapped spin liquid states have also 
been extensively studied [§4lfj|. 

Due to the difficulty of handling strong correlations, 
numerical simulations have been playing an important 
role on the study of exotic quantum spin states (lil - 

H EMU- 

Whether the spin-disordered quantum in- 
sulating states exist in the honeycomb lattice or not 
is currently under debating [3, [l6j]. A constrained 
path-integral quantum Monte-Carlo (QMC) simulation 
finds a gapless spin disordered phase in the square lat- 
tice with 7r-flux per plaquette [17] . Evidence of gapped 
spin liquid phases have also been found by the density- 
matrix-renormalization-group simulations of the frus- 
trated Heisenberg models in the Kagome lattice [l8| and 
in the square lattice with diagonal couplings jl9j |. 

The Fermi- Hubbard models with 27V components pos- 
sessing the SU(27V) or Sp(27V) symmetries are not only of 
academic interest now, but also have become the goal of 
experimental efforts in the ultra-cold atom physics [20j |. 
It was first proposed to use large-spin alkali and alkaline- 
earth atoms to realize the Sp(27V) and SU(27V) Hubbard 



models in Ref. [21] for the special case of 27V = 4 with the 
proof of a generic Sp(4) symmetry without fine-tuning. 
Currently, the SU(6) and SU(10) symmetric systems of 
173 Yb and 87 Sr atoms have been realized, respectively 

2^- 24 1 . In particular, the 173 Yb atoms have been loaded 
into optical lattices to realize the SU(6) Hubbard model, 
and the charge excitation gap has been observed [23j]. It 
has also been expected that Pomeranchuk cooling is effi- 
cient in the large- TV case to further cool the system down 
to the temperature scale of the AF exchanges [25[ . 

In this article, we investigate the magnetic properties 
of the half-filled SU(2TV) Hubbard models with 27V = 4 
and 6 by the sign-problem free determinant projector 
QMC method. For the SU(4) case, the ground state re- 
mains AF ordered as that of the SU(2) case although the 
residual spin moments are much weaker. Quantum spin 
fluctuations in the SU(6) case completely suppress the 
AF long-range order. Furthermore, crystalline orders in 
the spin singlet channel including the dimer and charge 
flux are also absent. The power-law decay of the AF cor- 
relations is consistent with a quantum disordered spin 
gapless phase. Different from the previous work in which 
the gapless spin disordered states are found in the weak- 
Mott insulating regime with two fermions per unit cell 

17] , our work shows such a state deeply inside the Mott- 
insulating phase with a large single-particle gap and an 
odd number of fermions per unit cell. 

The SU(27V) Fermi Hubbard model in the 2D square 
lattice at half-filling is defined as 



H = —t 



where t is scaled as 1 below; a represents spin indices 
running from 1 to 27V; denotes the summation over 
the nearest neighbors; rii is the particle number operator 
on site i defined as Hi = X)a=i C L C «"- Eq. [T]is invariant 
under the particle-hole transformation in bipartite lat- 
tices as Ci a — > (— ) l c| Q , and thus the average filling per 
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site (m) = N. Similarly to the case of SU(2), the SU(2JV) 
Hubbard model at half-filling in bipartite lattices is free 
of the sign problem for an arbitrary value of 2N. 

We use the determinant projector QMC method for 
fermions with the periodical boundary condition 26T 28j | . 
The simulated system sizes Lx L range from L = 4 to 16. 
Finite-size scaling is performed to extrapolate the ground 
state properties in the thermodynamic limit. The initial 
trial wavefunction is the ground state of the free part of 
Eq. [T] whose hopping integral is attached a small flux 
to break the degeneracy [151 ] ■ Such a Slater-determinant 
plane-wave state for the imaginary-time evolution is as- 
sumed to be non-orthogonal to the true ground state 
of the entire Hamiltonian. The second order Suzuki- 
Trotter decomposition is performed with the imaginary 
time-step At = 0.05. The convergence of the simulation 
results with respect to different values of At has been 
checked. The length of the imaginary-time evolution is 
(3 = 40. For the SU(2) case, the Hubbard-Stratonovich 
(HS) transformation is usually performed by using the 
discrete Ising spin fields [l[. However, the spin channel 
decomposition does not easily generalize to the SU (2N) 
case due to the increasing of spin components. Instead, 
we follow the approximate discrete HS decomposition in 
the density channel at the price of involving complex 
numbers [l|. The error of this approximation is at the 
order (At) 4 , smaller than that of the Suzuki- Trotter de- 
composition, thus is negligible. This method has the ad- 
vantage that the SU(2iV) symmetry is maintained explic- 
itly, and also it easily generalizes to large values of 2N. 

Let us fix the convention of the SU(2iV) generators. 
The Hilbert space on site i filled with r (1 < r < 2N) 
fermions forms the SU(2N) representation described by 
the single column Young pattern denoted as l r where r 
is the number of rows. For these ^-representations, the 
SU(2N) generators are defined as 



7=1 



4(i)c 7 (i). 



(2) 



Another standard definition is through the generalized 
Gell-mann matrices c^(i)A^C/3(i) with 1 < a < 47V 2 - 1 
and the normalization condition of tr[A°A''] = \5 ah ■ The 
definition in Eq. [2] has a simple commutation relation 
as [J"* 3 , J lS ] = 5 Pl J aS - 8 a sJ lf) - However, the price 
is that not all of the operators of Eq. [2] are indepen- 
dent, which satisfy the constraint J2 a J° a = ®- The 
quadratic Casimir operator is expressed as C2{2N) = 
\ J2up J a P{i)jP a (i)- For the V representation, its value 
is related to the filling number r through the Fierz iden- 
tity as C 2 (2N,r) = r{2N -r){2N + l)/(4/V). In the large 
U limit in which charge fluctuations are negligible, each 
site represents the self-conjugate representation 1^. The 
two-site equal time spin-spin correlation function is de- 



fined as 
C 



'j,su { 2N)(iJ) = c 2 (2N,N) ^l {jaP{l)J0a{j)) > (3) 

Q./3 



where C(2N,N) = N(2N + l)/4 is the Casimir for 1 N 
representation. Cj t su(2N)(hi) approaches 1 in the large 
?7-limit. The normalized spin structure factor at the AF 
wavevector Q is defined as 

s su(2N) (Q) = c J N ^ N) Y.\( ja ^Q) J " a (Q))' W 

where J a ^{Q) = ^^'^^^(j). The imaginary- 
time-displaced spin-spin correlations at wavevector Q are 
defined as 

Ssu(2N)(Q,r)=J2(J a HQ,r)J^(Q,0)}, (5) 

a/3 

which are used to extract spin gaps below. 
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FIG. 1: The half-filled SU(4) Hubbard model in the square 
lattice, (a) The appearance of the AF long-range-order from 
the finite size scaling of the spin structure factor at Q = (tt, 7r) 
for U — 6 and 8. Solid curves are quadratic fits of data. 
The inset shows a typical SU(4) AF configuration in which 
different colors represent different spin components, (b) The 
absence of the spin gap from the finite size scaling of A S (Q). 

The SU(4) case Below we present the study of quan- 
tum spin fluctuations starting with the SU(i) case in 
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the square lattice, in which we find long-range AF or- 
dering since intermediate values of U. The finite size 
scaling of the spin structure factor j^S su ^{Q) at the 
AF wavevector Q — (tt,tt) is plotted in Fig. H](a). For 
example, at U = 8, it extrapolates to a small but finite 
value of so = 0.025 as L — > oo, which indicates the ex- 
istence of the AF long-range Neel order. In comparison, 
for the SU(2) case at the same value of 7, the extrap- 
olated value of lim^^ -^S SU{2) {Q) w 0.118 0. This 
shows the enhancement of quantum spin fluctuations as 
2N increases. 

Let us bipartite the lattice into A and B sublattices. 
One typical classic SU(4) Neel configuration is that A- 
sites are filled with components 1 and 2, and -B-sites filled 
with components 3 and 4. SU(4) is a rank-3 Lie group, 
and thus its Cartan algebra has three commutable gen- 
erators defined as A'i^ = j^i( n i ~ "-2) ± ("-3 — "-4)], 
and A 3 = [{ni + n 2 ) — (n 3 + n 4 )]. Each site of 
the above 57(4) configuration is a singlet of A12, and 
with the eigenvalues of ±-^j for A 3 . The AF long- 
range-ordered states possess gapless Goldstone modes, 
and the Goldstone manifold is the 8-dimensional Grass- 
mann one 7(4)/[7(2) x (7(2)]. The spin excitations 
carry quantum numbers of A 12 .3 as (±-^,0, ±-75) and 



(a) 



(0, ±-^,±-^). To verify the absence of spin gap, we 
calculate the imaginary-time-displaced spin correlation 
function Ssu(i) {Qi T ) [21, [31|. The finite size spin-gap 



A s (Q, 1/L) is fitted from the slope of In Ssu(i) {Q, T ) v - s - 
t. The finite-size scaling is plotted in Fig. [T] (b) which 
shows the absence of spin gap in consistent with the AF 
ordering. 

The SU(6) case As 2N increases to 6, quantum spin 
fluctuations become even stronger. The QMC simulation 
of the spin structure factor at Q = (71", ir) is presented 
in Fig. [2] (a). The finite size scaling of the SU(6) AF 
structure factor shows the absence of long-range-order for 
all cases of 7 — 4, 8 and 12, which means that the ground 
state is quantum paramagnetic. We further calculate the 
spin gap value at Q = (7r, tt) from the imaginary-time- 
displaced 5(7(6) spin correlation function Ssu(6) {Q, r )> 
and plot the extracted spin gap values in Fig. [2] (b). The 
finite-size scaling shows the vanishing of spin gap in the 
SU(6) case for all the three values of U — 4, 8 and 12. 

The finite size scalings of the spin-structure factor and 
spin gap are consistent with a quantum spin paramag- 
netic phase. In order to understand its nature, we calcu- 
late the farthest two-point equal-time spin-spin correla- 
tion Cj i 5[/( 6 )(L/2, L/2), and plot it in Fig. [U (c) which 
shows an algebraic decay as Cj t su(6){L/2 7 L/2) L~ n . 
For the case of the largest value of 7 = 12, the decay 
power is r\ — 1.17. The fact that 77 > 1 is consistent 
with the absence of AF long-range order. If the AF were 
long-range ordered, the Goldstone modes have linear dis- 
persions which would give rise to the linear decay of the 
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FIG. 2: Spin correlations of the half-filled SU(6) Hubbard 
model, (a) The finite-size scaling of the spin structure fac- 
tor at Q = (tt, 7r) shows the absence of the AF long-range 
order at U — 4, 8 and 12. Solid curves are quadratic 
fits, (b) The finite-size scaling A S (Q) shows the absence of 
spin gap. (c) The scaling of the farthest point correlation 
C J}SU(6) (L/2,L/2) for U = 8 and 12. 



equal-time transverse spin-spin correlations. 

We also confirm this algebraic correlation from the be- 
havior of the AF structure factor. In addition to the 
quadratic fits as presented in Fig. [5] (a), we also fit it 
by using the form of Ssu(6){Q) = AL~ 2 + BL~ n where 
the L _2 -term represents the short-range correlation and 
L^-term represents the quasi-long-range correlation. As 
shown in the supplementary information, the fitted value 
of rj is consistent with that from the farthest two-point 
correlation. 
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FIG. 3: Spin singlet channel operators of the half-filled SU(6) 
Hubbard model, (a) The finite-size scaling of the columnar 
dimer structure factors at Q' = (n,0). (b) The finite-size 
scaling of the DDW structure factors at Q — (ji, n). 



We further check other possible ordering patterns in- 
volving two neighboring sites. At half-filling, total parti- 
cle number on a bond is 27V, which is sufficient to form 
an SU(2N) singlet to minimize the spin superexchange 
energy. We consider ordering patterns in the spin sin- 
glet channel with translational symmetry breaking. The 
bond dimer and current operators are defined as the real 
and imaginary parts of the hopping amplitudes between 
nearest neighbors as 



Dij — c i,a C j,a 



F, 



= £*(< 



h.c), (6) 



and d-density-wave (DDW) operators as DDW(i) = 
{—) l ^2jF{i,j) where fj — r\ — ±e x , and ±i y . In the 
large U limit, the Heisenberg term S a ^ {i)S^ a (j) is gen- 
erated from the second order virtual hopping process, 
thus Dij can be used as the dimer order parameter. The 
structure factor of at Q' = (it, 0) and that of DDW 
at Q = (tt, 71"), after divided by L 2 , and are plotted in Fig. 
[2]b) and c), respectively. They are fitted by a power-law 
(1/L) 2 , thus their correlations are short-ranged. 

Single-particle gap The single-particle gaps for the 
SU(4) and SU(6) Hubbard models are also calculated at 
half-filling through the onsite imaginary time-displaced 
Green's function G(0,r) = ^ Ei<*c|c(i, T)c+(i, 0)|* G ), 
where |>Fg) is the ground state. At long time displace- 
ment, G(0, r) — > e~ AcT where A c is the single-particle 



FIG. 4: Charge gaps of half-filled SU(27V) Hubbard models. 
A) Charge gaps with U = 8 at 27V = 2, 4 and 6. B) The 1/L 
scaling of the charge gap for the half-filled SU(6) model at 
[7=12. 



excitation gap, thus A c can be fitted from the slope of 
lnG(0,r) v.s. r. Let us consider the large [/-limit for 
an intuitive picture: in the Mott-insulating background, 
the energy of adding a particle is lowered from U by fur- 
ther virtual particle-hole excitations. In other words, the 
Mott-insulator is polarizable. As increasing 27V, the con- 
figuration numbers of the virtual particle-hole excitations 
increase, which enhances charge fluctuations and thus re- 
duces the single-particle gap. In Fig. H(a), A c s are plot- 
ted at a fixed [7 = 8 for 27V = 2,4 and 6, all of which 
are finite. The details of finite-size scaling for A c are 
presented in the supplementary material. For the SU(6) 
case, A c = 0.15 is rather small at U = 8. Nevertheless, 
A c increases to 1.26 at U — 12 at which the system is 
safely inside the Mott-insulating regime. The charge lo- 
calization length can be estimated as £ c w Vf / A c «3~4 
which is much smaller than the maximal sample size 
L = 16. 

In conclusion, we have studied the ground state quan- 
tum antiferromagnetism in a half-filled SU(27V) Hubbard 
model in square lattice. We find that at 27V = 4, a long- 
range AF order still survive. In contrast, at 27V = 6, 
we have checked the absence of AF, dimer, and charge 
flux orderings. Based on the absence of the spin gap 
and the power-law spin correlations, we conjecture that 
the ground state of half-filled SU(6) Hubbard model is a 
possible candidate for the U(l) algebraic spin liquid. Our 
SU(6) quantum disordered state exists at U = 12 with 
a large single-particle gap. Whether further increasing 
U will drive the system into the AF long-range-ordered 
phase is to be investigated. Furthermore, our system has 
three fermions per unit cell, which are different from two 
fermions per unit cell studied in Ref. fl5, 17[. 
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We present some technical details for the quantum 
Monte Carlo simulations for the half-filled SU(2iV) Hub- 
bard model in the square lattice. 



Method 

Hubbard- Stratonovich transformation: We follow the 
approximate discrete Hubbard-Stratonovich (HS) decom- 
position in the density channel at the price of involving 
complex numbers This advantage is that the SU(2Af) 
symmetry is maintained explicitly, and also it easily gen- 
eralizes to large values of 2N. The HS transformation for 
half-filled SU(2N) Hubbard model reads as: 



-A (m—N) _ 



E 

l=±l,±2 



^(lyniimm-N) +0 (A T 4) ; ( 7 ) 



where m = J2 a =i C L C "*> A = ^/AtII/2, and 
7(±D = l + f, 7(±2) = 1-^, 



(±1) = ±a/2(3-V6), t ? (±2)=± V /2(3 + V6). (8) 



The error of this approximation is at the order (At) 4 , 
smaller than that of the Suzuki- Trotter decomposition, 
thus is negligible. 

Projector Quantum Monte Carlo: At zero tempera- 
ture, the ground state (GS) wave function \^g) can 
be obtained by the projector quantum Monte Carlo 
(PQMC) method which projects a trial wave function 
\^t) in the following way: 



= lim 

/3->oo 



-0H 



* 7 



(9) 



where \^t) is required to be nonorthogonal to \^o). As 
(3 is large, this projection procedures can filter out states 
other than the GS. (3 plays the role as the projector pa- 
rameter. The expectation value in the zero temperature 
limit is defined as: 



(O) 



(*o|*o) 



(^ T | e -fgQ e -fg|^ r ) 

<*r|e-^|*T) ' 



(10) 



Similarly to the finite temperature scheme, the zero tem- 
perature problem can be formulated as the determi- 
nant QMC method with Suzuki- Trotter decompositions 
and auxiliary fields {/}, called projector quantum Monte 
Carlo (PQMC). 

Initial state: In the framework of the PQMC, the 
trial wavefunction is chosen as a ground state of the 
non-interacting Hamiltonian. We consider a free tight- 
binding Hamiltonian on the square lattice with a tiny 
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magnetic flux $ through the sample as 

H = -tJ2{c]c i e i ^^ dhA + h.c\, (11) 

where A = $e x /L is the vector potential and $o = hc/e. 
In Eq. [TTJ we neglect the spin index for simplicity. The 
purpose of introducing the flux is to break the ground 
state degeneracy at half-filling for Eq. [TT] Due to the 
SU(2iV) symmetry, the trial wavefunction for all the 2N 
components can be chosen the same, and thus the total 
trial wavefunction is the direct product of them. In our 
simulations, we take $/<I>o = 0.0001. 



The spin structure factor of the SU(6) Hubbard 
model 




4 6 8 10 12 14 16 L 



FIG. 5: Fitting using Eq. [12] for the spin structure factor 
Ssu(6){Q)/L 2 for SU(6) Hubbard model with U = 12. 

The critical spin liquid is characterized by the algebraic 
decay of the spin-spin correlations in the long distance. 
The spin structure factor S(Q)/L 2 , which also contains 
the short-range information of spin-spin correlations, can 
be fitted by the function: 

S(Q)/L 2 = A/L 2 + B/L v . (12) 

The first term and the second term come from the short- 
range and the long-range algebraic correlations, respec- 
tively. 

We use Eq. (JT2J} to fit the data of S(Q)/L 2 obtained 
by the QMC in Fig. 3 (a) (main text) as shown in Fig. 
[5j For values of a = 1.325, b = 0.531, and r\ = 1.11, the 
difference between the QMC and fitted results are within 
the error bars of the QMC simulation. The obtained 
exponent here r\ = 1.11 is consistent with that from the 
fit for the farthest correlations (r) = 1.17). 




FIG. 6: The imaginary-time-displaced on-site Green's func- 
tions for U = 12 with different sample sizes of L — 8, 12 and 
16. 
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FIG. 7: Finite scalings of the single-particle gap for the half- 
filled SU(2iV) Hubbard model with 2N = 2, 4 and 6 at U = 8. 
The extrapolated single-particle gap values are A c = 2.26, 
0.89, and 0.15, for 2N = 2,4 and 6, respectively. 



Derivations of the single-particle gaps from the 
time-displaced Green's functions 

We define the single-particle gap as the ground state 
energy change of adding a particle to the ground state of 
the iV-particle system as A c = E Q (N + 1) - E (N) - fi, 
where fj, accounts for the changing of the particle number. 
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The time-displaced Green's function is defined as 
G>(0,t) = (^(rHK) 

i 

= ^J2 e ~ T ^ +1 ' E °~^o\ci\^ +1 )\ 2 . 

i.n 

(13) 

Therefore, at large t, we have G > (f— 0, r) ~ e~ rAc to 
estimate the values of A c . 

We present the QMC simulation results for the 
imaginary-time-displaced Green's functions for U — 12 in 
Fig. [6] The slopes of lnG(0,r) v.s. r give rise the finite 



size single-particle gap presented in Fig. 4 (b) in the main 
text, in which the finite scaling shows that A c = 1.26. We 
also present the finite size scaling of the single-particle 
gap values presented in Fig. 4 (a) (main text) in Fig. [7] 
They are the single-particle gaps at U — 8 for half-filled 
SU(2iV) models with 2N = 2,4 and 6, which show the 
rapid decrease of gap values as increasing 27V. 



[1] F. F. Assaad, ArXiv: cond- mat/9806307 (1998). 
[2] F. F. Assaad and M. Imada, J. Phys. Soc. Jpn 65, 189 
(1996). 



